Baltimore_Food_Desert

Natalie Adams and Matthew Haraden

Story Memo Link: https://docs.google.com/document/d/1FNg8osWoG5mVxaWA0KETqtRD_43NwM_R2noyWTJPKB4/edit?usp=sharing

For our final project, we plan to focus on the prevalence of grocery stores in Baltimore City. Food deserts are an issue in cities, especially in low income and minority areas.

Looking into grocery store data in Baltimore can help us determine if these stores are more accessible to people of a higher income versus people of a lower income. It can also help us see if there are any particular grocery stores that appear more often in higher income communities. The same goes for majority white and non-white majority neighborhoods. This data can be used to find possible discrepancies for grocery store locations in Baltimore.

Questions: How many grocery stores are in low income neighborhoods? How many grocery stores are in high income neighborhoods? How many grocery stores are in majority non-white neighborhoods? How many grocery stores are in majority white neighborhoods? How does this compare to convenience stores?

Things to Consider: How to define high and low income How to define majority white and non-majority white

Grocery Datasets: Supermarkets, Small Grocery Stores and Convenience Stores in Maryland https://data-clf.hub.arcgis.com/datasets/cce1f648db4b41acafdd66e45d8ad529_595/explore?location=38.803090%2C-77.248650%2C8.00&showTable=true

Income Datasets: https://vital-signs-bniajfi.hub.arcgis.com/datasets/bniajfi::median-household-income-community-statistical-area/explore?location=39.284832%2C-76.620524%2C12.52&showTable=true

Race Datasets: https://vital-signs-bniajfi.hub.arcgis.com/datasets/bniajfi::total-population/explore?layer=0&location=39.284832%2C-76.620524%2C12.52&showTable=true https://vital-signs-bniajfi.hub.arcgis.com/datasets/bniajfi::percent-of-residents-white-caucasian-non-hispanic-community-statistical-area-1/explore?location=39.284832%2C-76.620524%2C12.52 https://vital-signs-bniajfi.hub.arcgis.com/datasets/bniajfi::percent-of-residents-black-african-american-non-hispanic-community-statistical-area/explore?location=39.284832%2C-76.620524%2C12.52 ​​https://vital-signs-bniajfi.hub.arcgis.com/datasets/bniajfi::percent-of-residents-hispanic-community-statistical-area-1/explore?location=39.284832%2C-76.620524%2C12.52 https://vital-signs-bniajfi.hub.arcgis.com/datasets/percent-of-residents-asian-non-hispanic-community-statistical-area-1/explore?location=39.284832%2C-76.620524%2C12.52&showTable=true https://vital-signs-bniajfi.hub.arcgis.com/datasets/bniajfi::percent-of-residents-all-other-races-hawaiian-pacific-islander-alaskan-native-american-other-race-non-hispanic-community-statistical-area-1/explore?location=39.284832%2C-76.620524%2C12.52 https://vital-signs-bniajfi.hub.arcgis.com/datasets/bniajfi::percent-of-residents-two-or-more-races-non-hispanic-community-statistical-area-1/explore?location=39.284832%2C-76.620524%2C12.52

#load libraries

# turn off scientific notation
options(scipen=999)
# load libraries
library(janitor)

Attaching package: ‘janitor’

The following objects are masked from ‘package:stats’:

    chisq.test, fisher.test

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ─────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.2      ✔ forcats 0.5.2 ── Conflicts ────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(sf)

Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(leaflet)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(magrittr)

Attaching package: ‘magrittr’

The following object is masked from ‘package:purrr’:

    set_names

The following object is masked from ‘package:tidyr’:

    extract

library(ggplot2)

#load data Grocery store data comes from the Johns Hopkins Center for a Livable Future Open Data Portal and median household income and race data comes from the Baltimore Neighborhood Indicators Alliance Vital Signs Open Data Portal (see links above).

# food stores in maryland
stores <- read.csv("food_stores.csv")
# median household income by baltimore neighborhood
income <- read.csv("median_household_income.csv")
# baltimore neighborhoods shapefile
neighborhoods <- st_read("income_shapefiles/median_household_income.shp")

Reading layer `median_household_income' from data source 
  `/Users/natalieadams/Desktop/Jour472_Fall_2022/Final_Project/income_shapefiles/median_household_income.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 55 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1393927 ymin: 557733.6 xmax: 1445503 ymax: 621406.8
Projected CRS: NAD83 / Maryland (ftUS)

# total population by baltimore neighborhood
total_pop1 <- read.csv("total_population.csv")
# white population by baltimore neighborhood
percent_white1 <- read.csv("percent_white.csv")
# black population by baltimore neighborhood
percent_black1 <- read.csv("percent_black.csv")
# hispanic population by baltimore neighborhood
percent_hispanic1 <- read.csv("percent_hispanic.csv")
# asian population by baltimore neighborhood
percent_asian1 <- read.csv("percent_asian.csv")
# other race population by baltimore neighborhood
percent_other1 <- read.csv("percent_other.csv")
# two or more races population by baltimore neighborhood
percent_more1 <- read.csv("percent_two_or_more.csv")
#baltimore population shapefile
population <- st_read("population_shapefiles/total_population.shp")

Reading layer `total_population' from data source 
  `/Users/natalieadams/Desktop/Jour472_Fall_2022/Final_Project/population_shapefiles/total_population.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 55 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1393927 ymin: 557733.6 xmax: 1445503 ymax: 621406.8
Projected CRS: NAD83 / Maryland (ftUS)

#clean race data Cleaned the data to select just the neighborhood and 2020 population and then joined all the races together into one table. We did this because we couldn’t locate a single dataframe organized by race and neighborhood, so we combined separate ones and checked to make sure the percentages added up to 100.

# total population by baltimore neighborhood clean
total_pop2 <- total_pop1 %>% 
  select(CSA2010, tpop20) %>% 
  clean_names() %>%
  rename(neighborhood = 1, total_pop_2020 = 2)
# white population by baltimore neighborhood clean
percent_white2 <- percent_white1 %>% 
  select(CSA2010, pwhite20) %>% 
  clean_names() %>%
  rename(neighborhood = 1, percent_white = 2)
# black population by baltimore neighborhood clean
percent_black2 <- percent_black1 %>% 
  select(CSA2010, paa20) %>% 
  clean_names() %>%
  rename(neighborhood = 1, percent_black = 2)
# hispanic population by baltimore neighborhood clean
percent_hispanic2 <- percent_hispanic1 %>% 
  select(CSA2010, phisp20) %>% 
  clean_names() %>%
  rename(neighborhood = 1, percent_hispanic = 2)
# asian population by baltimore neighborhood clean
percent_asian2 <- percent_asian1 %>% 
  select(CSA2010, pasi20) %>% 
  clean_names() %>%
  rename(neighborhood = 1, percent_asian = 2)
# other race population by baltimore neighborhood clean
percent_other2 <- percent_other1 %>% 
  select(CSA2010, ppac20) %>% 
  clean_names() %>%
  rename(neighborhood = 1, percent_other = 2)
# two or more races population by baltimore neighborhood clean
percent_more2 <- percent_more1 %>% 
  select(CSA2010, p2more20) %>% 
  clean_names() %>%
  rename(neighborhood = 1, percent_two_or_more = 2)
# total baltimore population by race and baltimore neighborhood clean
baltimore_population <- total_pop2 %>% 
  inner_join(percent_white2, by="neighborhood") %>% 
  inner_join(percent_black2, by="neighborhood") %>% 
  inner_join(percent_hispanic2, by="neighborhood") %>% 
  inner_join(percent_asian2, by="neighborhood") %>% 
  inner_join(percent_other2, by="neighborhood") %>% 
  inner_join(percent_more2, by="neighborhood")

#clean grocery store data Selected relevant data points, filtered it for just the Baltimore City stores and cleaned column names. Made separate dataframes for supermarkets, small grocery stores and convenience stores.

# baltimore stores
baltimore_stores <- stores %>% 
  select(Main_Categ, Company_Na, County, Address, ZIP_Code, X, Y) %>% 
  filter(County == "Baltimore City") %>% 
  clean_names() %>%
  rename(store_type = 1, name = 2) %>% 
  st_as_sf(coords = c("x", "y"), crs = 4269) %>% 
  filter(store_type != "Public Market")
# baltimore supermarkets
supermarkets <- read_csv("food_stores.csv") %>% 
  select(Main_Categ, Company_Na, County, Address, ZIP_Code, X, Y) %>% 
  filter(County == "Baltimore City") %>% 
  clean_names() %>%
  rename(store_type = 1, name = 2) %>%
  filter(store_type == "Supermarket") %>% 
  st_as_sf(coords = c("x", "y"), crs = 4269)

Rows: 4690 Columns: 12── Column specification ────────────────────────────────────────────
Delimiter: ","
chr (8): Main_Categ, Sub_Catego, Company_Na, Address, City, Stat...
dbl (4): X, Y, OBJECTID, ZIP_Code
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# baltimore small grocery stores
small_grocery <- read_csv("food_stores.csv") %>% 
  select(Main_Categ, Company_Na, County, Address, ZIP_Code, X, Y) %>% 
  filter(County == "Baltimore City") %>% 
  clean_names() %>%
  rename(store_type = 1, name = 2) %>%
  filter(store_type == "Small Grocery/Corner Store") %>% 
  st_as_sf(coords = c("x", "y"), crs = 4269)

Rows: 4690 Columns: 12── Column specification ────────────────────────────────────────────
Delimiter: ","
chr (8): Main_Categ, Sub_Catego, Company_Na, Address, City, Stat...
dbl (4): X, Y, OBJECTID, ZIP_Code
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# baltimore convenience stores
convenience_stores <- read_csv("food_stores.csv") %>% 
  select(Main_Categ, Company_Na, County, Address, ZIP_Code, X, Y) %>% 
  filter(County == "Baltimore City") %>% 
  clean_names() %>%
  rename(store_type = 1, name = 2) %>%
  filter(store_type == "Convenience Store") %>% 
  st_as_sf(coords = c("x", "y"), crs = 4269)

Rows: 4690 Columns: 12── Column specification ────────────────────────────────────────────
Delimiter: ","
chr (8): Main_Categ, Sub_Catego, Company_Na, Address, City, Stat...
dbl (4): X, Y, OBJECTID, ZIP_Code
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

#clean median household income data and neighborhoods shapefile Selected most recent median household income data and cleaned column names.

The median household income in Baltimore is $52,164, according to the U.S. Census Bureau. This is less than the national median household income of $65,000. We chose to use the Baltimore median household income instead of the national median household income because it is more representative of the city. We defined low income as a median household income below Baltimore’s $52,164 median household income and high income as above that number. There are 32 neighborhoods characterized as low income and 23 neighborhoods characterized as high income.

# Load and clean low income neighborhood shapefile
low_income <- st_read("income_shapefiles/median_household_income.shp") %>% 
  clean_names() %>%
  select(csa2010, mhhi20, geometry) %>% 
  rename(neighborhood = 1, median_income=2) %>%
  mutate(
    income_level = case_when(
        median_income < 52164 ~ "low",
        median_income >= 52164 ~ "high"
      )
  ) %>% 
  filter(income_level == "low") %>%
  st_transform(crs=4269)

Reading layer `median_household_income' from data source 
  `/Users/natalieadams/Desktop/Jour472_Fall_2022/Final_Project/income_shapefiles/median_household_income.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 55 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1393927 ymin: 557733.6 xmax: 1445503 ymax: 621406.8
Projected CRS: NAD83 / Maryland (ftUS)

# Load and clean high income neighborhood shapefile
high_income <- st_read("income_shapefiles/median_household_income.shp") %>% 
  clean_names() %>%
  select(csa2010, mhhi20, geometry) %>% 
  rename(neighborhood = 1, median_income=2) %>%
  mutate(
    income_level = case_when(
        median_income < 52164 ~ "low",
        median_income >= 52164 ~ "high"
      )
  ) %>% 
  filter(income_level == "high") %>%
  st_transform(crs=4269)

Reading layer `median_household_income' from data source 
  `/Users/natalieadams/Desktop/Jour472_Fall_2022/Final_Project/income_shapefiles/median_household_income.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 55 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1393927 ymin: 557733.6 xmax: 1445503 ymax: 621406.8
Projected CRS: NAD83 / Maryland (ftUS)

# table for leaflet
baltimore_neighborhoods_income <- neighborhoods %>% 
  clean_names() %>% 
  select(csa2010, mhhi20) %>% 
  rename(neighborhood = 1, median_income=2) %>% 
  mutate(
    income_level = case_when(
        median_income < 52164 ~ "low",
        median_income >= 52164 ~ "high"
      )
  )

#How does the number of supermarkets compare in low income and high income Baltimore neighborhoods? There are 25 supermarkets in low income neighborhoods and 24 supermarkets in high income neighborhoods. While low income neighborhoods have one more supermarket than high income neighborhoods, there are more low income neighborhoods (32) than high income neighborhoods (23), making the ratio of supermarkets to neighborhoods greater for high income neighborhoods. There are some low income neighborhoods without a supermarket at all, while there are enough supermarkets for every high income neighborhood, plus an additional one. This shows that high income residents have greater access to supermarkets and fresh food than low income residents.

# sort coordinates into low income neighborhoods
low_income_supermarkets <- supermarkets %>%
  st_join(low_income) %>% 
  filter(income_level == "low")
# sort coordinates into high income neighborhoods
high_income_supermarkets <- supermarkets %>%
  st_join(high_income) %>% 
  filter(income_level == "high")

#plot supermarkets in baltimore based on income

ggplot()+
geom_sf(data=neighborhoods)+
geom_sf(data=low_income_supermarkets %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=income_level))+
geom_sf(data=high_income_supermarkets %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=income_level))+
theme_minimal()+
labs(title = "Supermarkets in High and Low Income Baltimore Neighborhoods", 
       subtitle = "Johns Hopkins and Baltimore Neighborhood Indicators Alliance Data, 2022",
       caption = "Graphic by Natalie Adams")

#How does the number of small grocery and corner stores compare in low income and high income Baltimore neighborhoods? There are 289 small grocery/corner stores in low income neighborhoods and 115 small grocery/corner stores in high income neighborhoods. Low income neighborhoods have more than twice as many small grocery/corner stores than high income neighborhoods and have a greater ratio of small grocery/corner stores than high income neighborhoods. Small grocery/corner stores typically offer ethnic or cultural food options, such as Asian or South American. The prevalence of small grocery/corner stores in low income neighborhoods indicates a possible high non-white population in those neighborhoods. Small grocery/corner stores could also make up for the absence of access to supermarkets in these communities, although access to fresh food could still be limited.

# sort coordinates into low income neighborhoods
low_income_small_grocery <- small_grocery %>%
  st_join(low_income) %>% 
  filter(income_level == "low")
# sort coordinates into high income neighborhoods
high_income_small_grocery <- small_grocery %>%
  st_join(high_income) %>% 
  filter(income_level == "high")

#plot small grocery stores in baltimore based on income

ggplot()+
geom_sf(data=neighborhoods)+
geom_sf(data=low_income_small_grocery %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=income_level))+
geom_sf(data=high_income_small_grocery %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=income_level))+
theme_minimal()+
labs(title = "Small Grocery and Corner Stores in High and Low Income Baltimore Neighborhoods", 
       subtitle = "Johns Hopkins and Baltimore Neighborhood Indicators Alliance Data, 2022",
       caption = "Graphic by Natalie Adams")

#How does the number of convenience stores compare in low income and high income Baltimore neighborhoods? There are 166 convenience stores in low income neighborhoods and 116 convenience stores in high income neighborhoods. The average number of convenience stores in low and high income neighborhoods is virtually the same (5.1 in low income compared to 5.0 in high income). Convenience stores typically have “junk food” and snacks, like chips, soda and candy. Large chain convenience stores, such as Royal Farms and Wawa, have options to order hot food, like fried chicken and sandwiches. Some convenience stores have fresh fruit, such as apples and bananas, but overall, do not offer great access to fresh food. While there is an equal percentage of convenience stores in high income and low income neighborhoods, the lower percentage of supermarkets in low income neighborhoods could indicate a stronger reliance on convenience stores for buying food.

# sort coordinates into low income neighborhoods
low_income_convenience_stores <- convenience_stores %>%
  st_join(low_income) %>% 
  filter(income_level == "low")
# sort coordinates into high income neighborhoods
high_income_convenience_stores <- convenience_stores %>%
  st_join(high_income) %>% 
  filter(income_level == "high")

#plot convenience stores in baltimore based on income

ggplot()+
geom_sf(data=neighborhoods)+
geom_sf(data=low_income_convenience_stores %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=income_level))+
geom_sf(data=high_income_convenience_stores %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=income_level))+
theme_minimal()+
labs(title = "Convenience Stores in High and Low Income Baltimore Neighborhoods", 
       subtitle = "Johns Hopkins and Baltimore Neighborhood Indicators Alliance Data, 2022",
       caption = "Graphic by Natalie Adams")

#Income Leaflet Graph

pal <- colorFactor(
  palette = c('red', 'darkgreen', 'black'),
  domain = baltimore_stores$store_type)

pal2 <- colorFactor(
  palette = c('#fbd009', '#199bd6'),
  domain = baltimore_neighborhoods_income$income_level)

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data=low_income,
              color = ("#199bd6"),
              weight = 1.5,
              smoothFactor = 0.2,
              fillOpacity = 0.1,
              label = low_income$neighborhood) %>% 
  addPolygons(data=high_income,
              color = ("#fbd009"),
              weight = 1.5,
              smoothFactor = 0.2,
              fillOpacity = 0.1,
              label = high_income$neighborhood) %>%
  addCircles(data=baltimore_stores,
             weight = 2,
             color = ~pal(store_type),
             label= paste("Store Name:",baltimore_stores$name)) %>% 
    addLegend(
   position = "bottomleft",
    pal = pal,
    values = baltimore_stores$store_type,
    title = "Stores by Type", group="group_1") %>% 
  addLegend(
   position = "bottomright",
    pal = pal2,
    values = baltimore_neighborhoods_income$income_level,
    title = "Neighborhood Income", group="group_2") %>% 
  addLayersControl(overlayGroups = c("group_1","group_2"),
                   options = layersControlOptions(collapsed = FALSE))

Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'

#clean race data and neighborhoods shapefile Selected white data and cleaned column names.

We defined majority white neighborhoods as neighborhoods with a neighborhood with a 50% or grater white population and non majority white neighborhoods as a neighborhood with a less than 50% white population. There are 11 majority white neighborhoods and 44 non majority white neighborhoods.

# filter for majority white neighborhoods
majority_white <- st_read("income_shapefiles/median_household_income.shp") %>% 
  clean_names() %>%
  select(csa2010, geometry) %>% 
  rename(neighborhood = 1) %>% 
  inner_join(percent_white2, by="neighborhood") %>% 
  mutate(
    population_status = case_when(
        percent_white < 50 ~ "non_majority_white",
        percent_white >= 50 ~ "majority_white"
      )
  ) %>%
  filter(population_status == "majority_white") %>%
  st_transform(crs=4269)

Reading layer `median_household_income' from data source 
  `/Users/natalieadams/Desktop/Jour472_Fall_2022/Final_Project/income_shapefiles/median_household_income.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 55 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1393927 ymin: 557733.6 xmax: 1445503 ymax: 621406.8
Projected CRS: NAD83 / Maryland (ftUS)

# filter for non majority white neighborhoods
non_majority_white <- st_read("income_shapefiles/median_household_income.shp") %>% 
  clean_names() %>%
  select(csa2010, geometry) %>% 
  rename(neighborhood = 1) %>% 
  inner_join(percent_white2, by="neighborhood") %>% 
  mutate(
    population_status = case_when(
        percent_white < 50 ~ "non_majority_white",
        percent_white >= 50 ~ "majority_white"
      )
  ) %>%
  filter(population_status == "non_majority_white") %>%
  st_transform(crs=4269)

Reading layer `median_household_income' from data source 
  `/Users/natalieadams/Desktop/Jour472_Fall_2022/Final_Project/income_shapefiles/median_household_income.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 55 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1393927 ymin: 557733.6 xmax: 1445503 ymax: 621406.8
Projected CRS: NAD83 / Maryland (ftUS)

# table for leaflet
baltimore_neighborhoods_race <- baltimore_population %>% 
  select(neighborhood, percent_white) %>%
  mutate(
    population_status = case_when(
        percent_white < 50 ~ "non_majority_white",
        percent_white >= 50 ~ "majority_white"
      )
  )

#How does the number of supermarkets compare in majority white and non majority white Baltimore neighborhoods? There are 14 supermarkets in majority white neighborhoods and 35 supermarkets in non majority white neighborhoods. While non majority white neighborhoods have a higher number of supermarkets the ratio of supermarkets for majority white neighborhoods (14 supermarkets for 11 neighborhoods) is greater than the ratio of supermarkets for non majority white neighborhoods (35 supermarkets for 44 neighborhoods). This indicates that white residents have greater access to supermarkets and fresh food than non white residents.

# sort coordinates into majority white neighborhoods
majority_white_supermarkets <- supermarkets %>%
  st_join(majority_white) %>% 
  filter(population_status == "majority_white")
# sort coordinates into non majority white neighborhoods
non_majority_white_supermarkets <- supermarkets %>%
  st_join(non_majority_white) %>% 
  filter(population_status == "non_majority_white")

#plot supermarkets in baltimore based on race

ggplot()+
geom_sf(data=neighborhoods)+
geom_sf(data=majority_white_supermarkets %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=population_status))+
geom_sf(data=non_majority_white_supermarkets %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=population_status))+
theme_minimal()+
labs(title = "Supermarkets in Majority White and Non Majority White Baltimore Neighborhoods", 
       subtitle = "Johns Hopkins and Baltimore Neighborhood Indicators Alliance Data, 2022",
       caption = "Graphic by Natalie Adams")

#How does the number of small grocery and corner stores compare in majority white and non majority white Baltimore neighborhoods? There are 12 times as many small grocery/corner stores in non majority white neighborhoods (373) than majority white neighborhoods (31). Small grocery/corner stores typically offer ethnic or cultural food options, such as Asian or South American, so it makes sense for most of these stores to be in majority non white communities. It is unclear from this data whether the high number of small grocery/corner stores is due to preference of the neighborhood populations, a lack of access to other types of food stores (such as supermarkets), other factors or a combination of factors.

# sort coordinates into majority white neighborhoods
majority_white_small_grocery <- small_grocery %>%
  st_join(majority_white) %>% 
  filter(population_status == "majority_white")
# sort coordinates into non majority white neighborhoods
non_majority_white_small_grocery <- small_grocery %>%
  st_join(non_majority_white) %>% 
  filter(population_status == "non_majority_white")

#plot small grocery stores in baltimore based on race

ggplot()+
geom_sf(data=neighborhoods)+
geom_sf(data=majority_white_small_grocery %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=population_status))+
geom_sf(data=non_majority_white_small_grocery %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=population_status))+
theme_minimal()+
labs(title = "Small Grocery and Corner Stores in Majority White and Non Majority White Baltimore Neighborhoods", 
       subtitle = "Johns Hopkins and Baltimore Neighborhood Indicators Alliance Data, 2022",
       caption = "Graphic by Natalie Adams")

#How does the number of convenience stores compare in majority white and non majority white Baltimore neighborhoods? There are 49 convenience stores in majority white neighborhoods and 233 convenience stores in non majority white neighborhoods. The percentage of convenience stores in majority white and non majority white neighborhoods is nearly equal (4.5 in majority white compared to 5.3 in non majority white). Convenience stores typically have “junk food” and snacks, like chips, soda and candy. Large chain convenience stores, such as Royal Farms and Wawa, have options to order hot food, like fried chicken and sandwiches. Some convenience stores have fresh fruit, such as apples and bananas, but overall, do not offer great access to fresh food. While there is a nearly equal percentage of convenience stores in majority white and non majority white neighborhoods, the lower percentage of supermarkets in non majority white neighborhoods could indicate a stronger reliance on convenience stores for buying food.

# sort coordinates into majority white neighborhoods
majority_white_convenience_stores <- convenience_stores %>%
  st_join(majority_white) %>% 
  filter(population_status == "majority_white")
# sort coordinates into non majority white neighborhoods
non_majority_white_convenience_stores <- convenience_stores %>%
  st_join(non_majority_white) %>% 
  filter(population_status == "non_majority_white")

#plot convenience stores in baltimore based on race

ggplot()+
geom_sf(data=neighborhoods)+
geom_sf(data=majority_white_convenience_stores %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=population_status))+
geom_sf(data=non_majority_white_convenience_stores %>% st_as_sf(coords = c("x", "y"), crs = 4326), mapping=aes(color=population_status))+
theme_minimal()+
labs(title = "Convenience Stores in Majority White and Non Majority White Baltimore Neighborhoods", 
       subtitle = "Johns Hopkins and Baltimore Neighborhood Indicators Alliance Data, 2022",
       caption = "Graphic by Natalie Adams")

#leaflet race graph

pal1 <- colorFactor(
  palette = c('red', 'darkgreen', 'black'),
  domain = baltimore_stores$store_type)

pal2 <- colorFactor(
  palette = c('#fbd009', '#199bd6'),
  domain = baltimore_neighborhoods_race$population_status)

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data=majority_white,
              color = ("#fbd009"),
              weight = 1.5,
              smoothFactor = 0.2,
              fillOpacity = 0.1,
              label = majority_white$neighborhood) %>% 
  addPolygons(data=non_majority_white,
              color = ("#199bd6"),
              weight = 1.5,
              smoothFactor = 0.2,
              fillOpacity = 0.1,
              label = non_majority_white$neighborhood) %>%
  addCircles(data=baltimore_stores,
             weight = 2,
             color = ~pal(store_type),
             label= paste("Store Name:",baltimore_stores$name)) %>% 
    addLegend(
   position = "bottomleft",
    pal = pal1,
    values = baltimore_stores$store_type,
    title = "Stores by Type", group="group_1") %>% 
  addLegend(
   position = "bottomright",
    pal = pal2,
    values = baltimore_neighborhoods_race$population_status,
    title = "Neighborhood Race", group="group_2") %>% 
  addLayersControl(overlayGroups = c("group_1","group_2"),
                   options = layersControlOptions(collapsed = FALSE))

Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'